/**
 * @file random_util.h
 *
 * @date Sep 23, 2011
 * @author seth
 *
 * @brief This file is part of EBEX State Estimator, created for the EBEX project
 *
 * This software is copyright (C) 2011 Columbia University
 *
 * EBEX State Estimator is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * EBEX State Estimator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with EBEX State Estimator; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifndef RANDOM_UTIL_H
#define RANDOM_UTIL_H

#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include <mkl.h>

#include <ebex_math.h>
#include <ebex_util.h>
#include <ebex_lookup.h>

typedef enum
{
	RANDOM_MCG31 = VSL_BRNG_MCG31, /*!< A 31-bit multiplicative congruential generator. */
	RANDOM_R250 = VSL_BRNG_R250, /*!< A generalized feedback shift register generator. */
	RANDOM_MRG32K3A = VSL_BRNG_MRG32K3A, /*!< A combined multiple recursive generator with two components of order 3. */
	RANDOM_MCG59 = VSL_BRNG_MCG59, /*!< A 59-bit multiplicative congruential generator. */
	RANDOM_WH = VSL_BRNG_WH, /*!< A set of 273 Wichmann-Hill combined multiplicative congruential generators. */
	RANDOM_MT19937 = VSL_BRNG_MT19937, /*!< A Mersenne Twister pseudorandom number generator. */
	RANDOM_MT2203 = VSL_BRNG_MT2203, /*!< A set of 6024 Mersenne Twister pseudorandom number generators. */
	RANDOM_SFMT19937 = VSL_BRNG_SFMT19937, /*!< A SIMD-oriented Fast Mersenne Twister pseudorandom number generator. */
	RANDOM_SOBOL = VSL_BRNG_SOBOL, /*!< A 32-bit Gray code-based generator producing low-discrepancy sequences for dimensions 1 ≤ s ≤ 40; user-defined dimensions are also available. */
	RANDOM_NIEDERR = VSL_BRNG_NIEDERR, /*!< A 32-bit Gray code-based generator producing low-discrepancy sequences for dimensions 1 ≤ s ≤ 318; user-defined dimensions are also available. */
	RANDOM_IABSTRACT = VSL_BRNG_IABSTRACT, /*!< An abstract random number generator for integer arrays. */
	RANDOM_DABSTRACT = VSL_BRNG_DABSTRACT, /*!< An abstract random number generator for double precision floating-point arrays. */
	RANDOM_SABSTRACT = VSL_BRNG_SABSTRACT, /*!< An abstract random number generator for single precision floating-point arrays. */
} random_generator_t;

#define _random_dists(x,_)	\
	_(x, none)				\
	_(x, uniform)			\
	_(x, gaussian)			\
	_(x, gamma)				\
	_(x, dirichlet)			\
	_(x, beta)				\
	_(x, epanechnikov)		\
	_(x, pearsonii)			\
	_(x, spherical)			\
	_(x, binomial)			\
	_(x, weibull)

EBEX_LOOKUP_TABLE(random_dist, static);

typedef struct
{
	double			mean;
	double			stddev;
} gaussian_parameters_t;

typedef struct
{
	double			trials;
	double			probability;
} binomial_parameters_t;

typedef struct
{
	double			shape;
	double			scale;
} weibull_parameters_t;

typedef struct
{
	double			min;
	double			max;
} uniform_parameters_t;

typedef struct
{
	double			shape;
	double			scale;
} gamma_parameters_t;

typedef struct
{
	double			a;
	double			b;
} beta_parameters_t;

typedef struct
{
	double			shape;
} dirichlet_parameters_t;

typedef struct
{
	e_random_dist	type;

	union
	{
		double					raw_param[2];
		uniform_parameters_t	uniform;
		gaussian_parameters_t	gaussian;
		beta_parameters_t		beta;
		gamma_parameters_t		gamma;
		dirichlet_parameters_t	dirichlet;
		weibull_parameters_t	weibull;
		binomial_parameters_t	binomial;
	};

	double			offset;
	double			scale;
} random_sample_t;

typedef struct
{
	e_random_dist	type;

	union
	{
		uniform_parameters_t	uniform[3];
		gaussian_parameters_t	gaussian[3];
		beta_parameters_t		beta[3];
		gamma_parameters_t		gamma[3];
		dirichlet_parameters_t	dirichlet[3];
		weibull_parameters_t	weibull[3];
		binomial_parameters_t	binomial[3];
	};

	double			offset[3];
	double			scale[3];
} random_sample_3d_t;

bool ran_seed(const unsigned int seed);
void ran_init(void);
random_generator_t ran_get_generator(void);
bool ran_set_generator(random_generator_t m_generator);

double ran_sample_scale_offset(const random_sample_t *m_param, double m_scale, double m_offset);
void ran_sample_scale_offset_n(double *m_dest, const random_sample_t *m_param, double m_scale, double m_offset, size_t m_num, size_t m_stride);
void ran_sample_scale_offset_3d(const random_sample_3d_t *m_param, double *m_retvals, const double *m_scale, const double *m_offset);

int ran_int_uniform_range(const int min, const int max);
void ran_double_range_uniform_n(const double m_min, const double m_max, double *m_output, size_t m_num);
void ran_gaussian_n(const double m_mean, const double m_stddev, double *m_output, int m_num);
void ran_gaussian_md(const double *m_mean, const double *m_cov, double *m_output, int m_dim, int m_num);
void ran_gamma_n(const double m_shape, const double m_offset, const double m_scale, double *m_output, const int m_size);
void ran_int_permutation(const int m_size, int *m_inputarray, int *m_outputarray);
void ran_binomial_n(const int m_trials, const double m_probability, int *m_output, const int m_num);
void ran_weibull_n(const double m_shape, const double m_displacement, const double m_scale, double *m_output, int m_num);

double ran_sample_mean(random_sample_t *m_param);
double ran_sample_mode(random_sample_t *m_param);

static inline double ran_sample(const random_sample_t *m_param)
{
	return ran_sample_scale_offset(m_param, m_param->scale, m_param->offset);
}

static inline void ran_sample_n(double *m_dest, const random_sample_t *m_param, size_t m_num, size_t m_stride)
{
	ran_sample_scale_offset_n (m_dest, m_param, m_param->scale, m_param->offset, m_num, m_stride);
}

static inline void ran_sample_3d(const random_sample_3d_t *m_param, double *m_retvals)
{
	ran_sample_scale_offset_3d(m_param, m_retvals, m_param->scale, m_param->offset);
}

static inline int ran_int_uniform(void)
{
	return ran_int_uniform_range(INT_MIN, INT_MAX);
}

static inline double ran_gamma(const double m_shape, const double m_offset, const double m_scale)
{
	double retval;
	ran_gamma_n(m_shape, m_offset, m_scale, &retval, 1);
	return retval;
}
static inline void ran_dirichlet(const double m_alpha, const double *m_offset, double *m_output, const int m_dim)
{
	int i;
	double norm = 0.0;

	ran_gamma_n(m_alpha, 0.0, 1.0, m_output, m_dim);

	for (i = 0; i < m_dim; i++)
	{
		norm += m_output[i];
	}

	for (i = 0; i < m_dim; i++)
	{
		m_output[i] = m_output[i] / norm + m_offset[i];
	}
}

static inline void ran_dirichlet_n(const double m_alpha, double *m_output, const int m_dim, const int m_num)
{
	int i, j;
	int skip;
	double norm = 0.0;

	ran_gamma_n(m_alpha, 0.0, 1.0, m_output, m_dim * m_num);

	for (i = 0; i < m_num; i++)
	{
		norm = 0.0;
		skip = i * m_dim;
		for (j = 0; j < m_dim; j++)
			norm += m_output[skip + j];

		for (j = 0; j < m_dim; j++)
			m_output[skip + j] /= norm;
	}
}

static inline void ran_beta_n(const double a, const double b, double *m_output, const int m_size)
{
	double tempa[m_size];
	double tempb[m_size];

	ran_gamma_n(a, 0.0, 1.0, tempa, m_size);
	ran_gamma_n(b, 0.0, 1.0, tempb, m_size);

	vdAdd(m_size, tempa, tempb, tempb);
	vdDiv(m_size, tempa, tempb, m_output);
}

static inline double ran_beta(const double a, const double b)
{
	double retval;
	ran_beta_n(a, b, &retval, 1);
	return retval;
}

static inline double ran_epanechnikov(void)
{
	return ran_beta(2.0,2.0);
}

static inline void ran_epanechnikov_n(double *m_output, const int m_size)
{
	ran_beta_n(2.0, 2.0, m_output, m_size);
}

/**
 * Returns a uniformly distributed value on the unit sphere
 * @param m_output
 * @param m_scale
 * @param m_offset Can be NULL.  Otherwise vector
 * @param m_dim
 */
static inline void ran_unit_sphere_md(double *m_output, const double m_scale, const double *m_offset, const int m_dim)
{
	double norm;

	do
	{
		norm = 0.0;
		ran_gaussian_n(0.0, 1.0, m_output, m_dim);
		for (int i = 0; i < m_dim; i++)
		{
			norm += (m_output[i] * m_output[i]);
		}
	}while (norm > 1.0 || norm <= DBL_EPSILON);

	norm = m_scale / sqrt(norm);

	for (int i = 0; i < m_dim; i++)
	{
		m_output[i] = m_output[i] * norm;
	}
	if (m_offset)
	{
		for (int i = 0; i < m_dim; i++)
		{
			m_output[i] += m_offset[i];
		}
	}

}

static inline void ran_mvpearson_ii_md(double *m_output, const double *m_offset, const int m_dim)
{
	double scale_factor = sqrt(m_dim + 4) * sqrt(ran_beta(m_dim/2, 2));
	ran_unit_sphere_md(m_output, scale_factor, m_offset, m_dim);
}

static inline bool ran_bool(void)
{
	return (bool) (ran_int_uniform() <= 0);
}

static inline bool ran_bool_prob(const double m_prob)
{
	return (bool) ((unsigned)ran_int_uniform() <= (unsigned) (m_prob * (double) UINT_MAX));
}

/**
 * Returns an integer distributed uniformly on [0,max).  Max will never be returned.
 * Maximum #max value is INT_MAX
 * @param max upper bound on uniform distribution (non-inclusive)
 * @return random integer
 */
static inline int ran_uint_uniform(const int max)
{
	return ran_int_uniform_range(0, max);
}

static inline double ran_double_range_uniform(const double min, const double max)
{
	double retval;
	ran_double_range_uniform_n(min, max, &retval, 1);
	return retval;
}

static inline double ran_double(const double max)
{
	return ran_double_range_uniform(0.0, max);
}

static inline double ran_unit_uniform(void)
{
	return ran_double_range_uniform(0.0, 1.0);
}

static inline void ran_unit_uniform_n(double *m_output, size_t m_samples)
{
	ran_double_range_uniform_n(0.0, 1.0, m_output, m_samples);
}

static inline double ran_gaussian(const double m_mean, const double m_stddev)
{
	double retval;
	ran_gaussian_n(m_mean, m_stddev, &retval, 1);
	return retval;
}

static inline void ran_unit_gaussian_n(double *m_output, int m_num)
{
	ran_gaussian_n(0.0, 1.0, m_output, m_num);
}

#endif
